# include "scene.h"
# include <stdio.h>
# include <stdlib.h>
# include <math.h>

struct var_state
   {
    double ***q, qreq;
    float dxdef, dydef, dxreq, dyreq;
    unsigned short posx
   };


extern short cmtrx_1dim_interpol(double *x, unsigned short nxpt, double *xipl, unsigned short nxiplpt,\
 unsigned short norder, struct cmtrx *cm, double **tempmtrx);

static short interpol_2dim(struct var_state *var, struct cmtrx *cmx, struct cmtrx *cmy, short icell, short jcell);

#define NxCELL 5
#define NyCELL 4

#define nORDER 3

static double func(double x, double y);

int main()
 {
  double x[NxCELL+1], y[NyCELL+1], xmid[NxCELL], ymid[NyCELL];
  struct cmtrx cmx, cmy;
  double **v;
  short i,j,ecode;
  struct var_state var;

  v = (double **)malloc(sizeof(double *)*NxCELL);
  for(i = 0; i < NxCELL; i++)
    {
     v[i] = (double *)malloc(sizeof(double)*NyCELL);
//      for(j = 0; j < Ny-1; j++) v[i][j] = (double *)malloc(sizeof(double)*(Nz-1));
    }

  for(i = 0; i <= NxCELL; i++) x[i] = sqrt(i);
  for(j = 0; j <= NyCELL; j++) y[j] = sqrt(j);
//   for(k = 0; k < Nz; k++) z[k] = sqrt(k);
/*
  for(k = 0; k < Nz-1; k++)
     {
      zmid[k] = 0.5*(z[k]+z[k+1]);*/
      for(j = 0; j < NyCELL; j++)
         {
          ymid[j] = 0.5*(y[j]+y[j+1]);
          for(i = 0; i < NxCELL; i++)
             {
              xmid[i] = 0.5*(x[i]+x[i+1]);
              v[i][j] = func(x[i+1],ymid[j]);
//           printf("%E\t",v[i][j]);
             }
//        printf("\n");
         }
//      }


  cmx.respadi = 1;
  cmx.srcpadi = 1;
  cmx.respade = 0;
  cmx.srcpade = 0;
  ecode = cmtrx_1dim_interpol(&x[cmx.srcpadi], 1+NxCELL-cmx.srcpadi-cmx.srcpade,\
                              &xmid[cmx.respadi], NxCELL-cmx.respadi-cmx.respade, nORDER, &cmx, 0);

  printf("ecode = %d, nrows = %d, ncols = %d\n", ecode, cmx.nrows, cmx.ncols);
//     for(i = 0; i < cmx.nrows; i++) {for(j = 0; j <= cmx.iend[i]-cmx.istart[i]; j++) printf("%E\t",cmx.mtrx[i][j]); printf("\n");}


  cmy.respadi = 1;
  cmy.srcpadi = 0;
  cmy.respade = 1;
  cmy.srcpade = 0;
  ecode = cmtrx_1dim_interpol(&ymid[cmy.srcpadi], NyCELL-cmy.srcpadi-cmy.srcpade,\
                              &y[cmy.respadi], 1+NyCELL-cmy.respadi-cmy.respade, nORDER, &cmy, 0);

  printf("ecode = %d, nrows = %d, ncols = %d\n", ecode, cmy.nrows, cmy.ncols);
//    for(i = 0; i < cmy.nrows; i++) {for(j = 0; j <= cmy.iend[i]-cmy.istart[i]; j++) printf("%E\t",cmy.mtrx[i][j]); printf("\n");}

  var.q = v;
  var.dxdef = 0.5;
  var.dydef = 0;
  var.dxreq = 0;
  var.dyreq = 0.5;

  interpol_2dim(&var,&cmx,&cmy,4,2);

//   q = interpol_2dim(&cmy,&cmx,v,3,4,0.5,0);

  printf("Interpolated Value = %E, Real Value = %E\n", var.qreq, func(xmid[4], y[3]));

  for(i = 0; i < NxCELL; i++) free(v[i]); free(v);

  return 0;
 }


/* If the array **q is: q[nx][ny], the coefficient matrices should also be passed in the order
   of cmx, cmy. If this order is not maintained, the results will be errorneous.

   The icell, jcell speficied the cell-numbers at which the interpolation is to be done.
*/


static short interpol_2dim(struct var_state *var, struct cmtrx *cmx, struct cmtrx *cmy, short icell, short jcell)
   {
    double qx, qy;
    short istartx, iendx, istarty, iendy, i, j, indxi, indxj, ireq, jreq;

    ireq = icell+(short)(2.01*var->dxreq);// cmx->resflag);
    jreq = jcell+(short)(2.01*var->dyreq);

    qy = 0;
    istarty = cmy->istart[jreq-cmy->respadi];
    iendy = cmy->iend[jreq-cmy->respadi];

    for(indxj = istarty; indxj <= iendy; indxj++)
        {
         j = indxj+cmy->srcpadi-(int)(2.01*var->dydef);

         qx = 0;
         istartx = cmx->istart[ireq-cmx->respadi];
         iendx = cmx->iend[ireq-cmx->respadi];

         for(indxi = istartx; indxi <= iendx; indxi++)
            {
             i = indxi+cmx->srcpadi-(int)(2.01*var->dxdef);
             qx += cmx->mtrx[ireq-cmx->respadi][indxi-istartx]*var->q[i][j];
//  printf("i,j = %d,%d\n",i,j);
            }
//   printf("j = %d, qx = %E\n",j, qx);
          qy += cmy->mtrx[jreq-cmy->respadi][indxj-istarty]*qx;
         }

    var->qreq = qy;
    return 1;
   }


static double func(double x, double y)
  {
   return sqrt(x)*y;
  }
